/*************************************************************
** 	Incentive Effects of Recall Elections:  		        **
**	Evidence from Criminal Sentencing in California Courts 	**
** 	Analysis for submission.do								**
** 	Conduct all analysis for paper							**
** 	Latest revision: August 17, 2021				        **
************************************************************/

clear
set more off
set seed 10242004
set matsize 1000

cd "/Users/sanfordgordon/Dropbox/Recall Paper/Analysis/Program Files/JOP Replication"
use "recall_data_main.dta"

/********************************************************************************
*** Main Effect 		   													  ***
********************************************************************************/
/********************************************************************************
*** Graphical Evidence (Figures 1 in paper and C.1 in SI)                     ***
********************************************************************************/
/********************************************************************************
*** Unadjusted             										              ***
********************************************************************************/
local vrange = "0.2(0.1)0.5"
local fitmodel = "lfit"
local fitmodel2 = "lpoly"
local dv = "sentence_normed_truncated"
local ar = 0.67
/* Petition Announced */
preserve
keep if abs(time_a)<=45
fastxtile n_all = time_a, nq(50)
gen n = .
replace n = n_all
collapse (mean) `dv' (median) time_a, by(n)
gen weight = max(0,(45-abs(time_a))/45)
twoway 	   (`fitmodel2' `dv' time_a if time_a <= 0 , lcolor(gs8) lpattern(solid)) ///
	       (`fitmodel2' `dv' time_a if time_a >=0 , lcolor(gs8) lpattern(solid)) ///
		   (`fitmodel' `dv' time_a if time_a <= 0 [aweight=weight], lcolor(maroon) lpattern(dash)) ///
	       (`fitmodel' `dv' time_a if time_a >=0 [aweight=weight], lcolor(maroon) lpattern(dash)) ///
		   (scatter `dv' time_a, mcolor(black) msymbol(circle_hollow)) ///
		    ,  name(graph_announce_unadjusted_paper, replace) ///
			title("Petition Announcement, Unadjusted")  xtitle("Date") ///
			xlabel( 0 "Jun 6" -45 "Apr 22" 45 "Jul 21",angle(45) nogrid) ///
			ytitle("Normalized Sentence") ylabel(`vrange', format(%3.1f)  nogrid) ///
			xline(0, lcolor(black) lpattern(dash)) legend(off) scheme(plotplain) aspectratio(`ar')
restore

/* Recall Election*/
preserve
keep if abs(time_r)<=45
fastxtile n_all = time_r, nq(50)
gen n = .
replace n = n_all
collapse (mean) `dv' (median) time_r, by(n)
gen weight = max(0,(45-abs(time_r))/45)
twoway 	   (`fitmodel2' `dv' time_r if time_r <= 0 , lcolor(gs8) lpattern(solid)) ///
	       (`fitmodel2' `dv' time_r if time_r >=0 , lcolor(gs8) lpattern(solid)) ///
		   (`fitmodel' `dv' time_r if time_r <= 0 [aweight=weight], lcolor(maroon) lpattern(dash)) ///
	       (`fitmodel' `dv' time_r if time_r >=0 [aweight=weight], lcolor(maroon) lpattern(dash)) ///
		   (scatter `dv' time_r, mcolor(black) msymbol(circle_hollow)) ///
		    ,  name(graph_recall_unadjusted_paper, replace) ///
			title("Recall Election, Unadjusted")  xtitle("Date") ///
			xlabel( 0 "Jun 5" -45 "Apr 21" 45 "Jul 20",angle(45)  nogrid) ///
			ytitle("Normalized Sentence") ylabel(`vrange', format(%3.1f) nogrid) ///
			xline(0, lcolor(black) lpattern(dash)) legend(off) scheme(plotplain) aspectratio(`ar')
restore 

/* For Appendix: Signature Certification*/
preserve
keep if abs(time_s)<=45
fastxtile n_all = time_s, nq(50)
gen n = .
replace n = n_all
collapse (mean) `dv' (median) time_s, by(n)
gen weight = max(0,(45-abs(time_s))/45)
twoway 	   (`fitmodel2' `dv' time_s if time_s <= 0 , lcolor(gs8) lpattern(solid)) ///
	       (`fitmodel2' `dv' time_s if time_s >=0 , lcolor(gs8) lpattern(solid)) ///
		   (`fitmodel' `dv' time_s if time_s <= 0 [aweight=weight], lcolor(maroon) lpattern(dash)) ///
	       (`fitmodel' `dv' time_s if time_s >=0 [aweight=weight], lcolor(maroon) lpattern(dash)) ///
		   (scatter `dv' time_s, mcolor(black) msymbol(circle_hollow)) ///
		    ,  name(graph_sigs_unadjusted_paper, replace) ///
			title("Signature Certification, Unadjusted")  xtitle("Date") ///
			xlabel(0 "Jan 24" -45 "Dec 10" 45 "March 10",angle(45)  nogrid) ///
			ytitle("Normalized Sentence") ylabel(`vrange', format(%3.1f) nogrid) ///
			xline(0, lcolor(black) lpattern(dash)) legend(off) scheme(plotplain) aspectratio(`ar')
restore 

/********************************************************************************
***Residualized, using judge and statute fixed effects     			          ***
********************************************************************************/
/* Petition Announced */
preserve
drop if county == "Sacramento"
qui reg `dv' i.judge i.statute_code numberofcounts if abs(time_a)<=45,
qui predict ry if e(sample), resid
qui summ `dv' if e(sample), meanonly
qui replace ry = r(mean) + ry
qui reg time_a i.judge i.statute_code numberofcounts if abs(time_a)<=45
qui predict rx if e(sample), resid
qui summ time_a if e(sample), meanonly
qui replace rx = r(mean)+rx 

fastxtile n_all = rx, nq(50)
gen n = .
replace n = n_all

collapse (mean) ry (median) rx time_a, by(n)
gen weight = max(0,(45-abs(rx))/45)
twoway	   (`fitmodel2' ry rx if rx <= 0 , lcolor(gs8) lpattern(solid)) ///
	       (`fitmodel2' ry rx if rx  >=0 , lcolor(gs8) lpattern(solid)) ///
	       (`fitmodel' ry rx if rx <= 0 [aweight=weight], lcolor(maroon) lpattern(dash)) ///
	       (`fitmodel' ry rx if rx  >=0 [aweight=weight], lcolor(maroon) lpattern(dash)) ///
		   (scatter ry rx, mcolor(black) msymbol(circle_hollow)) ///
		    ,  name(graph_announce_adjusted_paper, replace) ///
			title("Petition Announcement, Adjusted")  xtitle("Date") ///
			xlabel( 0 "Jun 6" -45 "Apr 22" 45 "Jul 21",angle(45)  nogrid) ///
			ytitle("Normalized Sentence") ylabel(`vrange', format(%3.1f) nogrid) ///
			xline(0, lcolor(black) lpattern(dash)) legend(off) scheme(plotplain) aspectratio(`ar')			
restore

/* Recall Election*/
preserve
drop if county == "Sacramento"
qui reg `dv' i.judge i.statute_code numberofcounts if abs(time_r)<=45
qui predict ry if e(sample), resid
qui summ `dv' if e(sample), meanonly
qui replace ry = r(mean) + ry
qui reg time_r i.judge i.statute_code numberofcounts if abs(time_r)<=45
qui predict rx if e(sample), resid
qui summ time_r if e(sample), meanonly
qui replace rx = r(mean)+rx 
summ rx
fastxtile n_all = rx, nq(50)
gen n = .
replace n = n_all

collapse (mean) ry (median) rx time_r, by(n)
gen weight = max(0,(45-abs(rx))/45)
twoway 	   (`fitmodel2' ry rx if rx <= 0 , lcolor(gs8) lpattern(solid)) ///
	       (`fitmodel2' ry rx if rx  >=0 , lcolor(gs8) lpattern(solid)) ///
	       (`fitmodel' ry rx if rx <= 0 [aweight=weight], lcolor(maroon) lpattern(dash)) ///
	       (`fitmodel' ry rx if rx  >=0 [aweight=weight], lcolor(maroon) lpattern(dash)) ///
			(scatter ry rx, mcolor(black) msymbol(circle_hollow)) ///
		    ,  name(graph_recall_adjusted_paper, replace) ///
			title("Recall Election, Adjusted")  xtitle("Date") ///
			xlabel( 0 "Jun 5" -45 "Apr 21" 45 "Jul 20",angle(45) nogrid) ///
			ytitle("Normalized Sentence") ylabel(`vrange', format(%3.1f) nogrid) ///
			xline(0, lcolor(black) lpattern(dash)) legend(off) scheme(plotplain) aspectratio(`ar')
restore


/*Generate Figure 1 in the paper, and export as pdf */
gr combine graph_announce_unadjusted_paper graph_recall_unadjusted_paper ///
		   graph_announce_adjusted_paper graph_recall_adjusted_paper, ///
		   imargin(0 0 0 0) graphregion(margin(l=22 r=22)) name(all_rd_plots, replace)
graph export "combined_rd_plots.pdf", replace

/* For Appendix: Signature Certification*/
preserve
drop if county == "Sacramento"
qui reg `dv' i.judge i.statute_code numberofcounts if abs(time_s)<=45,
qui predict ry if e(sample), resid
qui summ `dv' if e(sample), meanonly
qui replace ry = r(mean) + ry
qui reg time_s i.judge i.statute_code numberofcounts if abs(time_s)<=45
qui predict rx if e(sample), resid
qui summ time_s if e(sample), meanonly
qui replace rx = r(mean)+rx 
summ rx
fastxtile n_all = rx, nq(50)
gen n = .
replace n = n_all

collapse (mean) ry (median) rx time_s, by(n)
gen weight = max(0,(45-abs(rx))/45)
twoway 	   (`fitmodel2' ry rx if rx <= 0 , lcolor(gs8) lpattern(solid)) ///
	       (`fitmodel2' ry rx if rx  >=0 , lcolor(gs8) lpattern(solid)) ///
	       (`fitmodel' ry rx if rx <= 0 [aweight=weight], lcolor(maroon) lpattern(dash)) ///
	       (`fitmodel' ry rx if rx  >=0 [aweight=weight], lcolor(maroon) lpattern(dash)) ///
			(scatter ry rx, mcolor(black) msymbol(circle_hollow)) ///
		    ,  name(graph_sigs_adjusted_paper, replace) ///
			title("Signature Certification, Adjusted")  xtitle("Date") ///
			xlabel( 0 "Jan 24" -45 "Dec 10" 45 "March 10",angle(45) nogrid) ///
			ytitle("Normalized Sentence") ylabel(`vrange', format(%3.1f) nogrid) ///
			xline(0, lcolor(black) lpattern(dash)) legend(off) scheme(plotplain) aspectratio(`ar')
restore

/*Generate Figure C.1 in the SI, and export as pdf*/
gr combine graph_sigs_unadjusted_paper graph_sigs_adjusted_paper, ///
imargin(0 0 0 0) graphregion(margin(l=22 r=22)) name(sig_cert_appendix, replace)

graph export "signatures_rd_plots.pdf", replace

/********************************************************************************
*** Table 1 (main results) 													  ***
********************************************************************************/
matrix table1 = J(6,4,.)
rdrobust sentence_normed_truncated time_a if abs(time_a)<180, kernel(tri) vce(cluster cluster_county_statute) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_rb_l)
matrix table1[1,1]=e(tau_bc)           /*Robust estimate*/
matrix table1[2,1]=e(se_tau_rb)        /*Robust standard error*/
matrix table1[3,1]= b_temp[1,1]        /*Left side intercept*/ 
matrix table1[4,1] = sqrt(v_temp[1,1]) /*Standard error of left side intercept*/
matrix table1[5,1]=e(h_l)              /* Bandwidth */
matrix table1[6,1]=e(N_h_l)+e(N_h_r)   /*Effective sample size */ 

rdrobust sentence_normed_truncated time_a if abs(time_a)<180, kernel(tri) vce(cluster cluster_judge_statute) ///
		  covs(numberofcounts s_* j_*) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_rb_l)
matrix table1[1,2]=e(tau_bc)           /*Robust estimate*/
matrix table1[2,2]=e(se_tau_rb)        /*Robust standard error*/
matrix table1[3,2]= b_temp[1,1]        /*Left side intercept*/ 
matrix table1[4,2] = sqrt(v_temp[1,1]) /*Standard error of left side intercept*/
matrix table1[5,2]=e(h_l)              /* Bandwidth */
matrix table1[6,2]=e(N_h_l)+e(N_h_r)   /*Effective sample size */ 

rdrobust sentence_normed_truncated time_r if abs(time_r)<180, kernel(tri) vce(cluster cluster_county_statute) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_rb_l)
matrix table1[1,3]=e(tau_bc)           /*Robust estimate*/
matrix table1[2,3]=e(se_tau_rb)        /*Robust standard error*/
matrix table1[3,3]= b_temp[1,1]        /*Left side intercept*/ 
matrix table1[4,3] = sqrt(v_temp[1,1]) /*Standard error of left side intercept*/
matrix table1[5,3]=e(h_l)              /* Bandwidth */
matrix table1[6,3]=e(N_h_l)+e(N_h_r)   /*Effective sample size */ 

rdrobust sentence_normed_truncated time_r if abs(time_r)<180, kernel(tri) vce(cluster cluster_judge_statute) ///
		  covs(numberofcounts s_* j_*) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_rb_l)
matrix table1[1,4]=e(tau_bc)           /*Robust estimate*/
matrix table1[2,4]=e(se_tau_rb)        /*Robust standard error*/
matrix table1[3,4]= b_temp[1,1]        /*Left side intercept*/ 
matrix table1[4,4] = sqrt(v_temp[1,1]) /*Standard error of left side intercept*/
matrix table1[5,4]=e(h_l)              /* Bandwidth */
matrix table1[6,4]=e(N_h_l)+e(N_h_r)   /*Effective sample size */ 

/*This is Table 1*/
matrix list table1


/********************************************************************************
*** Table 2a               													  ***
********************************************************************************/
matrix table2a = J(6,4,.)
rdrobust sentence_normed_to_atc time_a if abs(time_a)<180 & efiledate < mdy(6,6,2016), kernel(tri) vce(cluster cluster_county_statute) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_rb_l)
matrix table2a[1,1]=e(tau_bc)           /*Robust estimate*/
matrix table2a[2,1]=e(se_tau_rb)        /*Robust standard error*/
matrix table2a[3,1]= b_temp[1,1]        /*Left side intercept*/ 
matrix table2a[4,1] = sqrt(v_temp[1,1]) /*Standard error of left side intercept*/
matrix table2a[5,1]=e(h_l)              /* Bandwidth */
matrix table2a[6,1]=e(N_h_l)+e(N_h_r)   /*Effective sample size */ 

rdrobust sentence_normed_to_atc time_a if abs(time_a)<180  & efiledate < mdy(6,6,2016), kernel(tri) vce(cluster cluster_judge_statute) ///
		  covs(numberofcounts s_* j_*) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_rb_l)
matrix table2a[1,2]=e(tau_bc)           /*Robust estimate*/
matrix table2a[2,2]=e(se_tau_rb)        /*Robust standard error*/
matrix table2a[3,2]= b_temp[1,1]        /*Left side intercept*/ 
matrix table2a[4,2] = sqrt(v_temp[1,1]) /*Standard error of left side intercept*/
matrix table2a[5,2]=e(h_l)              /* Bandwidth */
matrix table2a[6,2]=e(N_h_l)+e(N_h_r)   /*Effective sample size */ 

rdrobust sentence_normed_to_atc time_r if abs(time_r)<180  & efiledate < mdy(6,5,2018), kernel(tri) vce(cluster cluster_county_statute) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_rb_l)
matrix table2a[1,3]=e(tau_bc)           /*Robust estimate*/
matrix table2a[2,3]=e(se_tau_rb)        /*Robust standard error*/
matrix table2a[3,3]= b_temp[1,1]        /*Left side intercept*/ 
matrix table2a[4,3] = sqrt(v_temp[1,1]) /*Standard error of left side intercept*/
matrix table2a[5,3]=e(h_l)              /* Bandwidth */
matrix table2a[6,3]=e(N_h_l)+e(N_h_r)   /*Effective sample size */ 

rdrobust sentence_normed_to_atc time_r if abs(time_r)<180  & efiledate < mdy(6,5,2018), kernel(tri) vce(cluster cluster_judge_statute) ///
		  covs(numberofcounts s_* j_*) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_rb_l)
matrix table2a[1,4]=e(tau_bc)           /*Robust estimate*/
matrix table2a[2,4]=e(se_tau_rb)        /*Robust standard error*/
matrix table2a[3,4]= b_temp[1,1]        /*Left side intercept*/ 
matrix table2a[4,4] = sqrt(v_temp[1,1]) /*Standard error of left side intercept*/
matrix table2a[5,4]=e(h_l)              /* Bandwidth */
matrix table2a[6,4]=e(N_h_l)+e(N_h_r)   /*Effective sample size */ 

/*This is Table 2a*/
matrix list table2a

/********************************************************************************
*** Table 2b               													  ***
********************************************************************************/
matrix table2b = J(6,4,.)
rdrobust charge_bargain time_a if abs(time_a)<180 & plea_guilty == 1 & efiledate < mdy(6,6,2016), kernel(tri) vce(cluster cluster_county_statute) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_rb_l)
matrix table2b[1,1]=e(tau_bc)           /*Robust estimate*/
matrix table2b[2,1]=e(se_tau_rb)        /*Robust standard error*/
matrix table2b[3,1]= b_temp[1,1]        /*Left side intercept*/ 
matrix table2b[4,1] = sqrt(v_temp[1,1]) /*Standard error of left side intercept*/
matrix table2b[5,1]=e(h_l)              /*Bandwidth */
matrix table2b[6,1]=e(N_h_l)+e(N_h_r)   /*Effective sample size */ 

rdrobust charge_bargain time_a if abs(time_a)<180 & plea_guilty == 1 & efiledate < mdy(6,6,2016), kernel(tri) vce(cluster cluster_judge_statute) ///
		  covs(numberofcounts s_* j_*) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_rb_l)
matrix table2b[1,2]=e(tau_bc)           /*Robust estimate*/
matrix table2b[2,2]=e(se_tau_rb)        /*Robust standard error*/
matrix table2b[3,2]= b_temp[1,1]        /*Left side intercept*/ 
matrix table2b[4,2] = sqrt(v_temp[1,1]) /*Standard error of left side intercept*/
matrix table2b[5,2]=e(h_l)              /*Bandwidth */
matrix table2b[6,2]=e(N_h_l)+e(N_h_r)   /*Effective sample size */ 

rdrobust charge_bargain time_r if abs(time_r)<180 & plea_guilty == 1 & efiledate < mdy(6,5,2018), kernel(tri) vce(cluster cluster_county_statute) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_rb_l)
matrix table2b[1,3]=e(tau_bc)           /*Robust estimate*/
matrix table2b[2,3]=e(se_tau_rb)        /*Robust standard error*/
matrix table2b[3,3]= b_temp[1,1]        /*Left side intercept*/ 
matrix table2b[4,3] = sqrt(v_temp[1,1]) /*Standard error of left side intercept*/
matrix table2b[5,3]=e(h_l)              /*Bandwidth */
matrix table2b[6,3]=e(N_h_l)+e(N_h_r)   /*Effective sample size */ 

rdrobust charge_bargain time_r if abs(time_r)<180 & plea_guilty == 1 & efiledate < mdy(6,5,2018), kernel(tri) vce(cluster cluster_judge_statute) ///
		  covs(numberofcounts s_* j_*) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_rb_l)
matrix table2b[1,4]=e(tau_bc)           /*Robust estimate*/
matrix table2b[2,4]=e(se_tau_rb)        /*Robust standard error*/
matrix table2b[3,4]= b_temp[1,1]        /*Left side intercept*/ 
matrix table2b[4,4] = sqrt(v_temp[1,1]) /*Standard error of left side intercept*/
matrix table2b[5,4]=e(h_l)              /*Bandwidth */
matrix table2b[6,4]=e(N_h_l)+e(N_h_r)   /*Effective sample size */ 

/*This is Table 2b*/
matrix list table2b


/********************************************************************************
*** Figure 2: Placebo Tests       										      ***
********************************************************************************/
/********************************************************************************
*** Top panels: No covariates 			      								  ***
********************************************************************************/
preserve
set more off
set matsize 2000
local startdate = mdy(1,1,2016)
local enddate = mdy(12,31,2016)
scalar length = 1+`enddate'-`startdate'
matrix results = J(length,6,.)

forval t = `startdate'(1)`enddate'{
	disp in yellow "Day: `t'"
	capture drop time_x dummy 
	gen time_x = edate - `t'
	gen dummy = edate >= `t'
	qui rdrobust sentence_normed_truncated time_x if abs(time_x)<180, ///
			  kernel(tri) vce(cluster cluster_county_statute) all
	local t0 = `t'-`startdate'+1
	matrix results[`t0',1]= e(tau_bc)
	matrix results[`t0',2]=2*normal(-abs(e(tau_bc)/e(se_tau_rb)))
	matrix results[`t0',3]= e(tau_bc) - 1.959964*e(se_tau_rb)
	matrix results[`t0',4]= e(tau_bc) + 1.959964*e(se_tau_rb)
	matrix results[`t0',5]=`t'
	matrix results[`t0',6] = e(b_l)
}

capture drop results*
svmat results	  
format results5 %tdMon_dd 
cumul results2, gen(cumul_unadjusted)

qui summ results2 if results5 == mdy(6,6,2016)
local arrow_x = r(mean)
twoway (line cumul_unadj cumul_unadj) ///
	   (scatter cumul_unadj results2 if results5 !=mdy(6,6,2016), mcolor(gs6) msymbol(circle_hollow)) ///
       (scatter cumul_unadj results2  if results5 == mdy(6,6,2016), mcolor(maroon) msymbol(diamond) msize(medium) ) ///
	   (pcarrowi 0.3 `arrow_x' 0.05 `arrow_x' (1), mcolor(maroon) lcolor(maroon) mfcolor(maroon) barbsize(1) msize(2)), text(0.34 0.025 "Actual", color(maroon)) ///
       yscale() name(gr_pvalues_unadj, replace) xtitle("Placebo p-value", size(small)) ///
	   ytitle("Empirical Cumulative Distribution", size(small)) legend(off) ///
	   ylab(,nogrid format(%3.1f)) xlab(,nogrid) ///
	   title("Placebo RD P-Values, Unadjusted",size(small))
capture drop n
sort cumul_unadj
qui count if results1 !=.
local n_temp = r(N)
gen n = _n
qui summ n if results5==mdy(6,6,2016)
local p_rank = r(mean)
disp (`n_temp'-`p_rank')/(`n_temp'-1)

local arrow_x =mdy(6,6,2016)
local limit = 45
twoway (rspike results3 results4 results5 if results5!=20611 & abs(results5-mdy(6,6,2016))<`limit',lcolor(gs11)) ///
		(scatter results1 results5 if results5!=20611 & abs(results5-mdy(6,6,2016))<`limit', mcolor(gs6) msymbol(circle_hollow)) ///
		(scatter results1 results5 if results5==20611 & abs(results5-mdy(6,6,2016))<`limit', mcolor(maroon) msymbol(diamond) ) ///
		(rspike results3 results4 results5 if results5==20611 & abs(results5-mdy(6,6,2016))<`limit',lcolor(maroon) lwidth(medthick)) ///
		(pcarrowi -0.1 `arrow_x' -0.02 `arrow_x' (1), mcolor(maroon) lcolor(maroon) mfcolor(maroon) barbsize(1) msize(2)), text(-.11 `arrow_x' "Actual", color(maroon)) ///
		xlabel(,angle(45) nogrid) yline(0,lcolor(black)) legend(off) name(gr_ts_unadj, replace) xtitle("Date",size(small)) ///
		ylab(,format(%3.1f) nogrid) ytitle("Estimate and 95% CI", size(small)) ///
		title("Placebo RD Estimates and 95% Confidence Intervals, Unadjusted",size(small)) 
restore

/*
*NOTE: Placebo Tests with Covariates take approximately 52 hours to run (tested on Intel(R) Core(TM) i5-4590S CPU processor with 16.0 GB RAM); 	
/********************************************************************************
*** Bottom Panels: With covariates 			      							  ***
********************************************************************************/

preserve
set more off
set matsize 2000
local startdate = mdy(1,1,2016)
local enddate = mdy(12,31,2016)
scalar length = 1+`enddate'-`startdate'
matrix results = J(length,6,.)

forval t = `startdate'(1)`enddate'{
	disp in yellow "Day: `t'"
	capture drop time_x dummy 
	gen time_x = edate - `t'
	gen dummy = edate >= `t'
	qui rdrobust sentence_normed_truncated time_x if abs(time_x)<180, ///
			  kernel(tri) vce(cluster cluster_county_statute) covs(numberofcounts s_* j_*) all
	local t0 = `t'-`startdate'+1
	matrix results[`t0',1]= e(tau_bc)
	matrix results[`t0',2]=2*normal(-abs(e(tau_bc)/e(se_tau_rb)))
	matrix results[`t0',3]= e(tau_bc) - 1.959964*e(se_tau_rb)
	matrix results[`t0',4]= e(tau_bc) + 1.959964*e(se_tau_rb)
	matrix results[`t0',5]=`t'
	matrix results[`t0',6] = e(b_l)
	
}
local ar = 0.67
capture drop results*
svmat results	  
format results5 %tdMon_dd 
capture drop cumul*
cumul results2, gen(cumul)
local ar = 0.67
qui summ results2 if results5 == mdy(6,6,2016)
local arrow_x = r(mean)
twoway (line cumul cumul) ///
	   (scatter cumul results2 if results5 !=mdy(6,6,2016), mcolor(gs6) msymbol(circle_hollow)) ///
       (scatter cumul results2  if results5 == mdy(6,6,2016), mcolor(maroon) msymbol(diamond) msize(medium) ) ///
	   (pcarrowi 0.3 `arrow_x' 0.05 `arrow_x' (1), mcolor(maroon) lcolor(maroon) mfcolor(red) barbsize(1) msize(2)), text(0.33 0.025 "Actual", color(maroon)) ///
       yscale() name(gr_pvalues_adj, replace) xtitle("Placebo p-value",size(medium)) ///
	   ytitle("Empirical Cumulative Distribution",size(medium)) legend(off) aspectratio(`ar') ///
	   ylab(,nogrid format(%3.1f)) xlab(,nogrid) ///
	   title("Placebo RD P-Values, Adjusted",size(medium))

capture drop results*
svmat results	  
format results5 %tdMon_dd 
local arrow_x = mdy(6,6,2016)
twoway (rspike results3 results4 results5 if results5!=20611,lcolor(gs11)) ///
		(scatter results1 results5 if results5!=20611, mcolor(gs6) msymbol(circle_hollow)) ///
		(scatter results1 results5 if results5==20611, mcolor(maroon) msymbol(diamond) ) ///
		(rspike results3 results4 results5 if results5==20611,lcolor(maroon) lwidth(medthick)) ///
		(pcarrowi 0.35 `arrow_x' 0.22 `arrow_x' (1), mcolor(maroon) lcolor(maroon) mfcolor(maroon) barbsize(1) msize(2)), text(0.38 `arrow_x' "Actual", color(maroon)) ///
		xlabel(,angle(45) nogrid) yline(0,lcolor(black)) legend(off) name(gr_ts_adj, replace) xtitle("Date",size(medium)) /// 
		ylab(-0.4(0.1)0.2,format(%3.1f) nogrid) ytitle("Estimate and 95% CI",size(medium)) aspectratio(`ar') ///
		title("Placebo RD Estimates and 95% Confidence Intervals, Adjusted",size(medium)) 
capture drop n
sort cumul
qui count if results1 !=.
local n_temp = r(N)
gen n = _n
qui summ n if results5==mdy(6,6,2016)
local p_rank = r(mean)
disp (`n_temp'-`p_rank')/(`n_temp'-1)		

gr combine gr_pvalues_unadj gr_ts_unadj gr_pvalues_adj gr_ts_adj, xsize(7.5) ysize(5)
graph export "placebo_tests_both.pdf", replace		
restore
*/


/********************************************************************************
*** Table 3 (Heterogeneous by Crime Type) 									  ***
********************************************************************************/

set more off
matrix table3 = J(6,6,.)
rdrobust sentence_normed_truncated time_a if abs(time_a)<180 & sex_crime==1, ///
kernel(tri) vce(cluster cluster_county_statute) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_cl_l)
matrix table3[1,1]=e(tau_bc) 			/*Robust estimate*/
matrix table3[2,1]=e(se_tau_rb) 		/*Robust standard error*/
matrix table3[3,1]= b_temp[1,1] 		/*Left side intercept*/ 
matrix table3[4,1] = sqrt(v_temp[1,1]) 	/*Standard error of left side intercept*/
matrix table3[5,1]=e(h_l) 				/*Bandwidth */
matrix table3[6,1]=e(N_h_l)+e(N_h_r) 	/*Effective sample size */ 

rdrobust sentence_normed_truncated time_a if abs(time_a)<180 & sex_crime==1, ///
kernel(tri) vce(cluster cluster_judge_statute) covs(numberofcounts s_* j_*) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_cl_l)
matrix table3[1,2]=e(tau_bc) 			/*Robust estimate*/
matrix table3[2,2]=e(se_tau_rb) 		/*Robust standard error*/
matrix table3[3,2]= b_temp[1,1] 		/*Left side intercept*/ 
matrix table3[4,2] = sqrt(v_temp[1,1]) 	/*Standard error of left side intercept*/
matrix table3[5,2]=e(h_l) 				/*Bandwidth */
matrix table3[6,2]=e(N_h_l)+e(N_h_r) 	/*Effective sample size */ 

rdrobust sentence_normed_truncated time_a if abs(time_a)<180 & violent_nonsex_667==1, ///
kernel(tri) vce(cluster cluster_county_statute) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_cl_l)
matrix table3[1,3]=e(tau_bc) 			/*Robust estimate*/
matrix table3[2,3]=e(se_tau_rb) 		/*Robust standard error*/
matrix table3[3,3]= b_temp[1,1] 		/*Left side intercept*/ 
matrix table3[4,3] = sqrt(v_temp[1,1]) 	/*Standard error of left side intercept*/
matrix table3[5,3]=e(h_l) 				/*Bandwidth */
matrix table3[6,3]=e(N_h_l)+e(N_h_r) 	/*Effective sample size */ 

rdrobust sentence_normed_truncated time_a if abs(time_a)<180 & violent_nonsex_667==1, ///
kernel(tri) vce(cluster cluster_judge_statute) covs(numberofcounts s_* j_*) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_cl_l)
matrix table3[1,4]=e(tau_bc) 			/*Robust estimate*/
matrix table3[2,4]=e(se_tau_rb) 		/*Robust standard error*/
matrix table3[3,4]= b_temp[1,1] 		/*Left side intercept*/ 
matrix table3[4,4] = sqrt(v_temp[1,1]) 	/*Standard error of left side intercept*/
matrix table3[5,4]=e(h_l) 				/*Bandwidth */
matrix table3[6,4]=e(N_h_l)+e(N_h_r) 	/*Effective sample size */ 

rdrobust sentence_normed_truncated time_a if abs(time_a)<180 & nonviolent_crime_667==1, ///
kernel(tri) vce(cluster cluster_county_statute) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_cl_l)
matrix table3[1,5]=e(tau_bc) 			/*Robust estimate*/
matrix table3[2,5]=e(se_tau_rb) 		/*Robust standard error*/
matrix table3[3,5]= b_temp[1,1] 		/*Left side intercept*/ 
matrix table3[4,5] = sqrt(v_temp[1,1]) 	/*Standard error of left side intercept*/
matrix table3[5,5]=e(h_l) 				/*Bandwidth */
matrix table3[6,5]=e(N_h_l)+e(N_h_r) 	/*Effective sample size */ 

rdrobust sentence_normed_truncated time_a if abs(time_a)<180 & nonviolent_crime_667==1, ///
kernel(tri) vce(cluster cluster_judge_statute) covs(numberofcounts s_* j_*) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_cl_l)
matrix table3[1,6]=e(tau_bc) 			/*Robust estimate*/
matrix table3[2,6]=e(se_tau_rb) 		/*Robust standard error*/
matrix table3[3,6]= b_temp[1,1] 		/*Left side intercept*/ 
matrix table3[4,6] = sqrt(v_temp[1,1]) 	/*Standard error of left side intercept*/
matrix table3[5,6]=e(h_l) 				/*Bandwidth */
matrix table3[6,6]=e(N_h_l)+e(N_h_r) 	/*Effective sample size */ 


/*This is Table 3*/
matrix list table3


/********************************************************************************
*** Table 4 (RD Estimates by Race of Defendant)							      ***
********************************************************************************/

set more off
matrix table4a = J(6,6,.)
matrix table4bc = J(12,2,.)

rdrobust sentence_normed_truncated time_a if abs(time_a)<180 & Black==1, ///
kernel(tri) vce(cluster cluster_county_statute) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_cl_l)
local black_rd_est_unadj =e(tau_bc)
local black_rd_se_unadj = (se_tau_rb)
local black_lsi_est_unadj = b_temp[1,1]
local black_lsi_se_unadj = sqrt(v_temp[1,1])
matrix table4a[1,1]=e(tau_bc) 			/*Robust estimate*/
matrix table4a[2,1]=e(se_tau_rb) 		/*Robust standard error*/
matrix table4a[3,1]= b_temp[1,1] 		/*Left side intercept*/ 
matrix table4a[4,1] = sqrt(v_temp[1,1]) /*Standard error of left side intercept*/
matrix table4a[5,1]=e(h_l) 				/*Bandwidth */
matrix table4a[6,1]=e(N_h_l)+e(N_h_r) 	/*Effective sample size */ 

rdrobust sentence_normed_truncated time_a if abs(time_a)<180 & Black==1, ///
kernel(tri) vce(cluster cluster_judge_statute) covs(numberofcounts s_* j_*) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_cl_l)
local black_rd_est_adj =e(tau_bc)
local black_rd_se_adj = (se_tau_rb)
local black_lsi_est_adj = b_temp[1,1]
local black_lsi_se_adj = sqrt(v_temp[1,1])
matrix table4a[1,2]=e(tau_bc) 			/*Robust estimate*/
matrix table4a[2,2]=e(se_tau_rb) 		/*Robust standard error*/
matrix table4a[3,2]= b_temp[1,1] 		/*Left side intercept*/ 
matrix table4a[4,2] = sqrt(v_temp[1,1]) /*Standard error of left side intercept*/
matrix table4a[5,2]=e(h_l) 				/*Bandwidth */
matrix table4a[6,2]=e(N_h_l)+e(N_h_r) 	/*Effective sample size */ 

rdrobust sentence_normed_truncated time_a if abs(time_a)<180 & Hispanic==1, ///
kernel(tri) vce(cluster cluster_county_statute) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_cl_l)
local hispanic_rd_est_unadj =e(tau_bc)
local hispanic_rd_se_unadj = (se_tau_rb)
local hispanic_lsi_est_unadj = b_temp[1,1]
local hispanic_lsi_se_unadj = sqrt(v_temp[1,1])
matrix table4a[1,3]=e(tau_bc) 			/*Robust estimate*/
matrix table4a[2,3]=e(se_tau_rb) 		/*Robust standard error*/
matrix table4a[3,3]= b_temp[1,1] 		/*Left side intercept*/ 
matrix table4a[4,3] = sqrt(v_temp[1,1]) /*Standard error of left side intercept*/
matrix table4a[5,3]=e(h_l) 				/*Bandwidth */
matrix table4a[6,3]=e(N_h_l)+e(N_h_r) 	/*Effective sample size */ 

rdrobust sentence_normed_truncated time_a if abs(time_a)<180 & Hispanic==1, ///
kernel(tri) vce(cluster cluster_judge_statute) covs(numberofcounts s_* j_*) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_cl_l)
local hispanic_rd_est_adj =e(tau_bc)
local hispanic_rd_se_adj = (se_tau_rb)
local hispanic_lsi_est_adj = b_temp[1,1]
local hispanic_lsi_se_adj = sqrt(v_temp[1,1])
matrix table4a[1,4]=e(tau_bc) 			/*Robust estimate*/
matrix table4a[2,4]=e(se_tau_rb) 		/*Robust standard error*/
matrix table4a[3,4]= b_temp[1,1] 		/*Left side intercept*/ 
matrix table4a[4,4] = sqrt(v_temp[1,1]) /*Standard error of left side intercept*/
matrix table4a[5,4]=e(h_l) 				/*Bandwidth */
matrix table4a[6,4]=e(N_h_l)+e(N_h_r) 	/*Effective sample size */ 

rdrobust sentence_normed_truncated time_a if abs(time_a)<180 & White==1, ///
kernel(tri) vce(cluster cluster_county_statute) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_cl_l)
local white_rd_est_unadj =e(tau_bc)
local white_rd_se_unadj = (se_tau_rb)
local white_lsi_est_unadj = b_temp[1,1]
local white_lsi_se_unadj = sqrt(v_temp[1,1])
matrix table4a[1,5]=e(tau_bc) 			/*Robust estimate*/
matrix table4a[2,5]=e(se_tau_rb) 		/*Robust standard error*/
matrix table4a[3,5]= b_temp[1,1] 		/*Left side intercept*/ 
matrix table4a[4,5] = sqrt(v_temp[1,1]) /*Standard error of left side intercept*/
matrix table4a[5,5]=e(h_l) 				/*Bandwidth */
matrix table4a[6,5]=e(N_h_l)+e(N_h_r) 	/*Effective sample size */ 

rdrobust sentence_normed_truncated time_a if abs(time_a)<180 & White==1, ///
kernel(tri) vce(cluster cluster_judge_statute) covs(numberofcounts s_* j_*) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_cl_l)
local white_rd_est_adj =e(tau_bc)
local white_rd_se_adj = (se_tau_rb)
local white_lsi_est_adj = b_temp[1,1]
local white_lsi_se_adj = sqrt(v_temp[1,1])
matrix table4a[1,6]=e(tau_bc) 			/*Robust estimate*/
matrix table4a[2,6]=e(se_tau_rb) 		/*Robust standard error*/
matrix table4a[3,6]= b_temp[1,1] 		/*Left side intercept*/ 
matrix table4a[4,6] = sqrt(v_temp[1,1]) /*Standard error of left side intercept*/
matrix table4a[5,6]=e(h_l) 				/*Bandwidth */
matrix table4a[6,6]=e(N_h_l)+e(N_h_r) 	/*Effective sample size */ 



matrix table4bc[1,1]=abs(`white_rd_est_unadj'-`black_rd_est_unadj')
matrix table4bc[2,1]=sqrt((`white_rd_se_unadj')^2+(`black_rd_se_unadj')^2)
matrix table4bc[3,1]=abs(`hispanic_rd_est_unadj'-`black_rd_est_unadj')
matrix table4bc[4,1]=sqrt((`hispanic_rd_se_unadj')^2+(`black_rd_se_unadj')^2)
matrix table4bc[5,1]=abs(`white_rd_est_unadj'-`hispanic_rd_est_unadj')
matrix table4bc[6,1]=sqrt((`white_rd_se_unadj')^2+(`hispanic_rd_se_unadj')^2)
matrix table4bc[7,1]=abs(`white_lsi_est_unadj'-`black_lsi_est_unadj')
matrix table4bc[8,1]=sqrt((`white_lsi_se_unadj')^2+(`black_lsi_se_unadj')^2)
matrix table4bc[9,1]=abs(`hispanic_lsi_est_unadj'-`black_lsi_est_unadj')
matrix table4bc[10,1]=sqrt((`hispanic_lsi_se_unadj')^2+(`black_lsi_se_unadj')^2)
matrix table4bc[11,1]=abs(`white_lsi_est_unadj'-`hispanic_lsi_est_unadj')
matrix table4bc[12,1]=sqrt((`white_lsi_se_unadj')^2+(`hispanic_lsi_se_unadj')^2)

matrix table4bc[1,2]=abs(`white_rd_est_adj'-`black_rd_est_adj')
matrix table4bc[2,2]=sqrt((`white_rd_se_adj')^2+(`black_rd_se_adj')^2)
matrix table4bc[3,2]=abs(`hispanic_rd_est_adj'-`black_rd_est_adj')
matrix table4bc[4,2]=sqrt((`hispanic_rd_se_adj')^2+(`black_rd_se_adj')^2)
matrix table4bc[5,2]=abs(`white_rd_est_adj'-`hispanic_rd_est_adj')
matrix table4bc[6,2]=sqrt((`white_rd_se_adj')^2+(`hispanic_rd_se_adj')^2)
matrix table4bc[7,2]=abs(`white_lsi_est_adj'-`black_lsi_est_adj')
matrix table4bc[8,2]=sqrt((`white_lsi_se_adj')^2+(`black_lsi_se_adj')^2)
matrix table4bc[9,2]=abs(`hispanic_lsi_est_adj'-`black_lsi_est_adj')
matrix table4bc[10,2]=sqrt((`hispanic_lsi_se_adj')^2+(`black_lsi_se_adj')^2)
matrix table4bc[11,2]=abs(`white_lsi_est_adj'-`hispanic_lsi_est_adj')
matrix table4bc[12,2]=sqrt((`white_lsi_se_adj')^2+(`hispanic_lsi_se_adj')^2)

/*This is Table 4a */
matrix list table4a

/*This is Table 4b and 4c*/
matrix list table4bc

/********************************************************************************
*** Table 5 (Indirect Racial Burden Hypothesis)							      ***
********************************************************************************/

matrix table5 = J(7,6,.)

reg max_possible_sentence Black Hispanic  if time_a < 0  & (White==1|Black==1|Hispanic==1), cluster(county)
matrix b_temp = e(b) 
matrix v_temp = e(V)
matrix table5[1,1] = b_temp[1,1]
matrix table5[2,1] = sqrt(v_temp[1,1])
matrix table5[3,1] = b_temp[1,2]
matrix table5[4,1] = sqrt(v_temp[2,2])
matrix table5[5,1] = b_temp[1,3]
matrix table5[6,1] = sqrt(v_temp[3,3])
matrix table5[7,1]=e(N)

xtreg max_possible_sentence Black Hispanic if time_a < 0  & (White==1|Black==1|Hispanic==1), fe i(judge) cluster(judge)
matrix b_temp = e(b) 
matrix v_temp = e(V)
matrix table5[1,2] = b_temp[1,1]
matrix table5[2,2] = sqrt(v_temp[1,1])
matrix table5[3,2] = b_temp[1,2]
matrix table5[4,2] = sqrt(v_temp[2,2])
matrix table5[5,2] = b_temp[1,3]
matrix table5[6,2] = sqrt(v_temp[3,3])
matrix table5[7,2]=e(N)

reg max_possible_sentence Black Hispanic  if time_a >= 0  & (White==1|Black==1|Hispanic==1), cluster(county)
matrix b_temp = e(b) 
matrix v_temp = e(V)
matrix table5[1,3] = b_temp[1,1]
matrix table5[2,3] = sqrt(v_temp[1,1])
matrix table5[3,3] = b_temp[1,2]
matrix table5[4,3] = sqrt(v_temp[2,2])
matrix table5[5,3] = b_temp[1,3]
matrix table5[6,3] = sqrt(v_temp[3,3])
matrix table5[7,3]=e(N)

xtreg max_possible_sentence Black Hispanic if time_a >= 0  & (White==1|Black==1|Hispanic==1), fe i(judge) cluster(judge)
matrix b_temp = e(b) 
matrix v_temp = e(V)
matrix table5[1,4] = b_temp[1,1]
matrix table5[2,4] = sqrt(v_temp[1,1])
matrix table5[3,4] = b_temp[1,2]
matrix table5[4,4] = sqrt(v_temp[2,2])
matrix table5[5,4] = b_temp[1,3]
matrix table5[6,4] = sqrt(v_temp[3,3])
matrix table5[7,4]=e(N)

reg max_possible_sentence Black Hispanic if (White==1|Black==1|Hispanic==1), cluster(county)
matrix b_temp = e(b) 
matrix v_temp = e(V)
matrix table5[1,5] = b_temp[1,1]
matrix table5[2,5] = sqrt(v_temp[1,1])
matrix table5[3,5] = b_temp[1,2]
matrix table5[4,5] = sqrt(v_temp[2,2])
matrix table5[5,5] = b_temp[1,3]
matrix table5[6,5] = sqrt(v_temp[3,3])
matrix table5[7,5]=e(N)

xtreg max_possible_sentence Black Hispanic if (White==1|Black==1|Hispanic==1) , fe i(judge) cluster(judge)
matrix b_temp = e(b) 
matrix v_temp = e(V)
matrix table5[1,6] = b_temp[1,1]
matrix table5[2,6] = sqrt(v_temp[1,1])
matrix table5[3,6] = b_temp[1,2]
matrix table5[4,6] = sqrt(v_temp[2,2])
matrix table5[5,6] = b_temp[1,3]
matrix table5[6,6] = sqrt(v_temp[3,3])
matrix table5[7,6]=e(N)

/*This is Table 5*/

matrix list table5

/********************************************************************************
*** Figure 3 (Effect of Petition Announcement on Raw Sentence Length by Race) ***
********************************************************************************/

matrix results = J(6,4,.)
local condlist = "White Black Hispanic"
forval i = 1(1)3 {
	local cond: word `i' of `condlist'
	rdrobust log10_sentence_days time_a if abs(time_a)<180 & `cond' == 1, ///
	kernel(tri) vce(cluster cluster_county_statute) all
	matrix vl = e(V_rb_l)
	matrix vr = e(V_rb_r)
	matrix results[2*`i'-1,1] =  e(tau_bc_l)
	matrix results[2*`i'-1,2] =  e(tau_bc_l) - 1.959964*sqrt(vl[1,1])
	matrix results[2*`i'-1,3] =  e(tau_bc_l) + 1.959964*sqrt(vl[1,1])
	matrix results[2*`i'-1,4]=0
	matrix results[2*`i',1] =  e(tau_bc_r)
	matrix results[2*`i',2] = e(tau_bc_r) - 1.959964*sqrt(vr[1,1])
	matrix results[2*`i',3] = e(tau_bc_r) + 1.959964*sqrt(vr[1,1])
	matrix results[2*`i',4]=1
		}
capture drop results*	
svmat results
capture drop r
gen r = "White" in 1/2
replace r = "Black" in 3/4
replace r = "Hispanic" in 5/6
twoway (rspike results2 results3 results4) (scatter results1 results4) (line results1 results4), ///
	   by(r, col(3) legend(off) note("")) xtitle("") ytitle("Sentence Length (Days)") ///
	   xlabel(0 "Pre" 1 "Post") ylabel(2 "100" 2.39794 "250" 2.69897 "500" 2.8750613 "750")
graph export "sentence_length_by_race.pdf", replace 	
drop results*

/********************************************************************************
*** Figure 4 (Aggregate Effects Estimates) 									  ***
********************************************************************************/
capture drop agg_results*
matrix agg_results = J(5,3,.)
/********************************************************************************
*** Assume LATE is ATE                                                    ***
********************************************************************************/
/********************************************************************************
*** (1a). Unadjusted                                                          ***
********************************************************************************/
summ max_possible if announce == 1 & abs(time_a)<45
scalar N_window = r(N)
scalar mean_window = r(mean)
rdrobust sentence_days time_a if abs(time_a)<180, kernel(tri) vce(cluster cluster_county_statute) all
/*Note -- these results are on the normalized scale -- adjusted later in program*/
matrix agg_results[1,1] = (e(tau_bc)-invnorm(.975)*e(se_tau_rb))*N_window/365
matrix agg_results[1,2] = e(tau_bc)*N_window/365
matrix agg_results[1,3] = (e(tau_bc)+invnorm(.975)*e(se_tau_rb))*N_window/365

/********************************************************************************
*** (1b). Adjusted                                                            ***
********************************************************************************/
rdrobust sentence_days time_a if abs(time_a)<180, kernel(tri) /// 
		 vce(cluster cluster_judge_statute) covs(numberofcounts s_* j_*)
/*Note -- these results are on the normalized scale -- adjusted later in program*/
matrix agg_results[2,1] = (e(tau_bc)-invnorm(.975)*e(se_tau_rb))*N_window/365
matrix agg_results[2,2] = e(tau_bc)*N_window/365
matrix agg_results[2,3] = (e(tau_bc)+invnorm(.975)*e(se_tau_rb))*N_window/365

/********************************************************************************
*** 2. Parametric prediction                                                 ***
********************************************************************************/
/**NOTE: This uses a uniform kernel to get better estimates on the slopes to the left and right of the cutpoint.
This has minimal effect on the RD estimates. 
 ***/
capture drop s_hat*
capture drop diff
set matsize 1000
matrix param_results = J(500,1,.)
forval i=1(1)500 {
	preserve
	qui keep if abs(time_a)<45
	bsample _N
	qui reg sentence_normed_truncated announce##c.time_a i.judge i.statute_code if abs(time_a)<45
	qui predict s_hat1 if e(sample) & time_a>=0 & time_a<45
	/* Counterfactual predictions; assuming time trend of zero*/
	qui replace announce = 0
	qui replace time_a =0
	qui predict s_hat0 if e(sample) & time_a>=0 & time_a<45
	qui gen diff = max_possible*(s_hat1-s_hat0)
	qui summ diff
	scalar effect = r(mean)*r(N)/365
	disp in yellow effect
	matrix param_results[`i',1]=effect
	restore
}
capture drop param_results*
svmat param_results
centile param_results, centile(2.5 50 97.5)
matrix agg_results[3,1] = r(c_1)
matrix agg_results[3,2] = r(c_2)
matrix agg_results[3,3] = r(c_3)

/********************************************************************************
*** CIA-Based Estimators                                                      ***
********************************************************************************/
/*Determine sample based on positivity restrictions*/
preserve
keep if time_a > -40 & time_a < 45
egen positivity_judge = mean(announce), by(judge)
egen positivity_statute = mean(announce), by(statute_code)
summ positivity*
keep if positivity_judge > 0 & positivity_judge < 1 & positivity_statute > 0 & positivity_statute<1
count
/*Test Conditional independence assumptions */
qui reg sentence_normed_truncated i.judge i.statute_code numberofcounts time_a if announce == 0 
lincom time_a
qui reg sentence_normed_truncated  i.judge i.statute_code numberofcounts time_a if announce == 1 
lincom time_a
restore

/********************************************************************************
*** (3). Linear Reweighting Estimator										  ***
********************************************************************************/
capture program drop linreweight
program  linreweight, rclass
	version 14.2
	syntax varlist(numeric fv), t(varname numeric) m(varname numeric)
	tokenize `varlist'
	qui reg `1' `2' `3' if `t' == 0
	tempvar s_hat0
	qui _predict `s_hat0' if `t' == 1, xb
	qui reg `1' `2' `3' if `t' == 1
	qui tempvar s_hat1
	qui _predict `s_hat1' if  `t' == 1, xb
	qui tempvar diff
	qui gen `diff' = (`s_hat1'-`s_hat0')*`m'
	qui summ `diff'
	return scalar d = r(mean)*r(N)/365
	scalar d1 = r(mean)*r(N)/365
	disp in yellow d1
end


matrix lrwe = J(500,1,.)
forval i=1(1)500 {
	preserve
	qui keep if time_a>-40 & time_a<45
	qui egen positivity_judge = mean(announce), by(judge)
	qui egen positivity_statute = mean(announce), by(statute_code)
	qui keep if positivity_judge > 0 & positivity_judge < 1 & positivity_statute > 0 & positivity_statute<1
	bsample _N
	gen one = 1
	linreweight sentence_normed_truncated i.judge i.statute_code numberofcounts, t(announce) m(max_possible)
	matrix lrwe[`i',1]=r(d)
	restore
	}	
capture drop lrwe
svmat lrwe
centile lrwe, centile(2.5 50 97.5)
matrix agg_results[4,1] = r(c_1)
matrix agg_results[4,2] = r(c_2)
matrix agg_results[4,3] = r(c_3)


/********************************************************************************
*** (4) Propensity Score Estimator                                            ***
********************************************************************************/
matrix pse_results = J(500,1,.)
forval i=1(1)500 {
	preserve
	qui keep if time_a>-40 & time_a<45
	qui egen positivity_judge = mean(announce), by(judge)
	qui egen positivity_statute = mean(announce), by(statute_code)
	qui keep if positivity_judge > 0 & positivity_judge < 1 & positivity_statute > 0 & positivity_statute<1
	qui logit announce i.judge i.statute_code numberofcounts
	predict lambda, p
	bsample _N
	qui count if announce == 1
	local n_t1 = r(N)
	qui summ announce if e(sample)
	local tmean = r(mean)
	qui gen y_cf = sentence_normed_truncated*max_possible*(announce-lambda)/(`tmean'*(1-lambda)) 
	qui summ y_cf
	matrix pse_results[`i',1]=r(mean)*`n_t1'/365
	drop y_cf
	restore
	}
capture drop pse_results
svmat pse_results
centile pse_results, centile(2.5 50 97.5)
matrix agg_results[5,1] = r(c_1)
matrix agg_results[5,2] = r(c_2)
matrix agg_results[5,3] = r(c_3)


preserve
qui keep if time_a>-40 & time_a<45
tab announce
qui egen positivity_judge = mean(announce), by(judge)
qui egen positivity_statute = mean(announce), by(statute_code)
qui keep if positivity_judge > 0 & positivity_judge < 1 & positivity_statute > 0 & positivity_statute<1
tab announce
restore
	
/*This is Figure 4 */
capture drop agg_results*	
svmat agg_results
capture drop n
gen n = _n in 1/5
replace n = 0.5 in 6
replace n = 5.5 in 7
twoway  (rspike agg_results1 agg_results3 n , lcolor(gs11) lwidth(medthick)) (scatter agg_results2 n, msymbol(circle_hollow) mcolor(gs6) ), ///
	xtitle("") legend(off) xlabel(1 `" "(1a)" "Unadjusted LATE" "Extrapolation" "' 2 `" "(1b)" "Adjusted LATE" "Extrapolation" "' /// 
	3 `" "(2)" "Fully" "Parametric" "' 4 `" "(3)" "Linear" "Reweighting" "' 5 `" "(4)" "Propensity" "Score" "', ///
	labsize(medsmall)) ytitle(Additional Years Incarceration)	

graph export "aggregate_effects.pdf", replace	
l agg_results* in 1/6



/********************************************************************************
*** Table A1. Descriptive Statistics										  ***
********************************************************************************/
gen male = arrest_gender_def=="Male"
replace male =. if arrest_gender_def=="NA"|mi(arrest_gender_def)

capture drop tablea1*
matrix tablea1 = J(13,6,.)

local varlist = "sentence_days sentence_normed_to_statute sentence_normed_truncated sex_crime violent_nonsex_667 nonviolent_crime_667 Black Hispanic White male arrest_age_def"
forvalues i = 1/11 {
	local var : word `i' of `varlist'
	qui sum `var'
	matrix tablea1[`i',1]=`r(mean)'
	matrix tablea1[`i',2]=`r(sd)'
	qui sum `var' if abs(time_a)<45
	matrix tablea1[`i',3]=`r(mean)'
	matrix tablea1[`i',4]=`r(sd)'
	qui sum `var' if abs(time_r)<45
	matrix tablea1[`i',5]=`r(mean)'
	matrix tablea1[`i',6]=`r(sd)'
	}

*number of cases
gen nvals=1
sum nvals
matrix tablea1[12,1]=`r(sum)'
matrix tablea1[12,2]=`r(sum)'
sum nvals if abs(time_a)<45
matrix tablea1[12,3]=`r(sum)'
matrix tablea1[12,4]=`r(sum)'
sum nvals if abs(time_r)<45
matrix tablea1[12,5]=`r(sum)'
matrix tablea1[12,6]=`r(sum)'

*number of defendants
drop nvals
by defendant, sort: gen nvals = _n == 1 
replace nvals =. if mi(defendant)
sum nvals
matrix tablea1[13,1]=`r(sum)'
matrix tablea1[13,2]=`r(sum)'
sum nvals if abs(time_a)<45
matrix tablea1[13,3]=`r(sum)'
matrix tablea1[13,4]=`r(sum)'
sum nvals if abs(time_r)<45
matrix tablea1[13,5]=`r(sum)'
matrix tablea1[13,6]=`r(sum)'

/*This is Table A1*/
matrix list tablea1

/********************************************************************************
*** Figure A.1. Media Coverage        									      ***
********************************************************************************/
/* See Figure A.1.R and associated data                                        */


/********************************************************************************
*** Table B.1. Replication Using Uncensored Normalized Sentences    		  ***
********************************************************************************/
matrix tableb1 = J(6,4,.)
rdrobust sentence_normed_to_statute time_a if abs(time_a)<180, kernel(tri) vce(cluster cluster_county_statute) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_rb_l)
matrix tableb1[1,1]=e(tau_bc)           /*Robust estimate*/
matrix tableb1[2,1]=e(se_tau_rb)        /*Robust standard error*/
matrix tableb1[3,1]= b_temp[1,1]        /*Left side intercept*/ 
matrix tableb1[4,1] = sqrt(v_temp[1,1]) /*Standard error of left side intercept*/
matrix tableb1[5,1]=e(h_l)              /* Bandwidth */
matrix tableb1[6,1]=e(N_h_l)+e(N_h_r)   /*Effective sample size */ 

rdrobust sentence_normed_to_statute time_a if abs(time_a)<180, kernel(tri) vce(cluster cluster_judge_statute) ///
		  covs(numberofcounts s_* j_*) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_rb_l)
matrix tableb1[1,2]=e(tau_bc)           /*Robust estimate*/
matrix tableb1[2,2]=e(se_tau_rb)        /*Robust standard error*/
matrix tableb1[3,2]= b_temp[1,1]        /*Left side intercept*/ 
matrix tableb1[4,2] = sqrt(v_temp[1,1]) /*Standard error of left side intercept*/
matrix tableb1[5,2]=e(h_l)              /* Bandwidth */
matrix tableb1[6,2]=e(N_h_l)+e(N_h_r)   /*Effective sample size */ 

rdrobust sentence_normed_to_statute time_r if abs(time_r)<180, kernel(tri) vce(cluster cluster_county_statute) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_rb_l)
matrix tableb1[1,3]=e(tau_bc)           /*Robust estimate*/
matrix tableb1[2,3]=e(se_tau_rb)        /*Robust standard error*/
matrix tableb1[3,3]= b_temp[1,1]        /*Left side intercept*/ 
matrix tableb1[4,3] = sqrt(v_temp[1,1]) /*Standard error of left side intercept*/
matrix tableb1[5,3]=e(h_l)              /* Bandwidth */
matrix tableb1[6,3]=e(N_h_l)+e(N_h_r)   /*Effective sample size */ 

rdrobust sentence_normed_to_statute time_r if abs(time_r)<180, kernel(tri) vce(cluster cluster_judge_statute) ///
		  covs(numberofcounts s_* j_*) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_rb_l)
matrix tableb1[1,4]=e(tau_bc)           /*Robust estimate*/
matrix tableb1[2,4]=e(se_tau_rb)        /*Robust standard error*/
matrix tableb1[3,4]= b_temp[1,1]        /*Left side intercept*/ 
matrix tableb1[4,4] = sqrt(v_temp[1,1]) /*Standard error of left side intercept*/
matrix tableb1[5,4]=e(h_l)              /* Bandwidth */
matrix tableb1[6,4]=e(N_h_l)+e(N_h_r)   /*Effective sample size */ 

/*This is Table B.1*/
matrix list tableb1



/********************************************************************************
*** Table B.2. Replication Using Non-Normalized Sentence Length	    		  ***
********************************************************************************/
matrix tableb2 = J(6,4,.)
rdrobust sentence_days time_a if abs(time_a)<180, kernel(tri) vce(cluster cluster_county_statute) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_rb_l)
matrix tableb2[1,1]=e(tau_bc)           /*Robust estimate*/
matrix tableb2[2,1]=e(se_tau_rb)        /*Robust standard error*/
matrix tableb2[3,1]= b_temp[1,1]        /*Left side intercept*/ 
matrix tableb2[4,1] = sqrt(v_temp[1,1]) /*Standard error of left side intercept*/
matrix tableb2[5,1]=e(h_l)              /* Bandwidth */
matrix tableb2[6,1]=e(N_h_l)+e(N_h_r)   /*Effective sample size */ 

rdrobust sentence_days time_a if abs(time_a)<180, kernel(tri) vce(cluster cluster_judge_statute) ///
		  covs(numberofcounts s_* j_*) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_rb_l)
matrix tableb2[1,2]=e(tau_bc)           /*Robust estimate*/
matrix tableb2[2,2]=e(se_tau_rb)        /*Robust standard error*/
matrix tableb2[3,2]= b_temp[1,1]        /*Left side intercept*/ 
matrix tableb2[4,2] = sqrt(v_temp[1,1]) /*Standard error of left side intercept*/
matrix tableb2[5,2]=e(h_l)              /* Bandwidth */
matrix tableb2[6,2]=e(N_h_l)+e(N_h_r)   /*Effective sample size */ 

rdrobust sentence_days time_r if abs(time_r)<180, kernel(tri) vce(cluster cluster_county_statute) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_rb_l)
matrix tableb2[1,3]=e(tau_bc)           /*Robust estimate*/
matrix tableb2[2,3]=e(se_tau_rb)        /*Robust standard error*/
matrix tableb2[3,3]= b_temp[1,1]        /*Left side intercept*/ 
matrix tableb2[4,3] = sqrt(v_temp[1,1]) /*Standard error of left side intercept*/
matrix tableb2[5,3]=e(h_l)              /* Bandwidth */
matrix tableb2[6,3]=e(N_h_l)+e(N_h_r)   /*Effective sample size */ 

rdrobust sentence_days time_r if abs(time_r)<180, kernel(tri) vce(cluster cluster_judge_statute) ///
		  covs(numberofcounts s_* j_*) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_rb_l)
matrix tableb2[1,4]=e(tau_bc)           /*Robust estimate*/
matrix tableb2[2,4]=e(se_tau_rb)        /*Robust standard error*/
matrix tableb2[3,4]= b_temp[1,1]        /*Left side intercept*/ 
matrix tableb2[4,4] = sqrt(v_temp[1,1]) /*Standard error of left side intercept*/
matrix tableb2[5,4]=e(h_l)              /* Bandwidth */
matrix tableb2[6,4]=e(N_h_l)+e(N_h_r)   /*Effective sample size */ 

/*This is Table B.2*/
matrix list tableb2

/********************************************************************************
*** Table B.3. Replication using One Count Only 							  ***
********************************************************************************/
matrix tableb3 = J(6,4,.)
rdrobust sentence_normed_truncated time_a if abs(time_a)<180&numberofcounts==1, kernel(tri) vce(cluster cluster_county_statute) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_rb_l)
matrix tableb3[1,1]=e(tau_bc)           /*Robust estimate*/
matrix tableb3[2,1]=e(se_tau_rb)        /*Robust standard error*/
matrix tableb3[3,1]= b_temp[1,1]        /*Left side intercept*/ 
matrix tableb3[4,1] = sqrt(v_temp[1,1]) /*Standard error of left side intercept*/
matrix tableb3[5,1]=e(h_l)              /* Bandwidth */
matrix tableb3[6,1]=e(N_h_l)+e(N_h_r)   /*Effective sample size */ 

rdrobust sentence_normed_truncated time_a if abs(time_a)<180&numberofcounts==1, kernel(tri) vce(cluster cluster_judge_statute) ///
		  covs(s_* j_*) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_rb_l)
matrix tableb3[1,2]=e(tau_bc)           /*Robust estimate*/
matrix tableb3[2,2]=e(se_tau_rb)        /*Robust standard error*/
matrix tableb3[3,2]= b_temp[1,1]        /*Left side intercept*/ 
matrix tableb3[4,2] = sqrt(v_temp[1,1]) /*Standard error of left side intercept*/
matrix tableb3[5,2]=e(h_l)              /* Bandwidth */
matrix tableb3[6,2]=e(N_h_l)+e(N_h_r)   /*Effective sample size */ 

rdrobust sentence_normed_truncated time_r if abs(time_r)<180&numberofcounts==1, kernel(tri) vce(cluster cluster_county_statute) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_rb_l)
matrix tableb3[1,3]=e(tau_bc)           /*Robust estimate*/
matrix tableb3[2,3]=e(se_tau_rb)        /*Robust standard error*/
matrix tableb3[3,3]= b_temp[1,1]        /*Left side intercept*/ 
matrix tableb3[4,3] = sqrt(v_temp[1,1]) /*Standard error of left side intercept*/
matrix tableb3[5,3]=e(h_l)              /* Bandwidth */
matrix tableb3[6,3]=e(N_h_l)+e(N_h_r)   /*Effective sample size */ 

rdrobust sentence_normed_truncated time_r if abs(time_r)<180&numberofcounts==1, kernel(tri) vce(cluster cluster_judge_statute) ///
		  covs(s_* j_*) all
matrix b_temp = e(beta_p_l) 
matrix v_temp = e(V_rb_l)
matrix tableb3[1,4]=e(tau_bc)           /*Robust estimate*/
matrix tableb3[2,4]=e(se_tau_rb)        /*Robust standard error*/
matrix tableb3[3,4]= b_temp[1,1]        /*Left side intercept*/ 
matrix tableb3[4,4] = sqrt(v_temp[1,1]) /*Standard error of left side intercept*/
matrix tableb3[5,4]=e(h_l)              /* Bandwidth */
matrix tableb3[6,4]=e(N_h_l)+e(N_h_r)   /*Effective sample size */ 

/*This is Table B.3*/
matrix list tableb3

/********************************************************************************
*** Figure B.1. Balance on Charge Indicators (Petition Only) 				  ***
********************************************************************************/

set more off
matrix results = J(382,7,.)

forval t = 1(1)382{
	capture drop x dv
	disp in yellow "Charge: `t'"
	*count number of charges per day
	qui bys time_a: egen x=count(s_`t') if s_`t'==1
	qui replace x = 0 if mi(x)
	by time_a: egen dv = max(x)
	qui sum max_possible if s_`t'==1, d
	if missing("`r(max)'")==1{
		local max_poss = .
	}
	else if missing("`r(max)'")==0{
		local max_poss = `r(max)'
	}
	
	*estimate RD
	preserve
	keep time_a dv announce
	qui duplicates drop _all, force
	qui rdbwselect dv time_a if abs(time_a)<180, kernel(tri)
	scalar bw = e(h_mserd)
	capture drop weight
	qui gen weight = max(0,(bw-abs(time_a))/bw)
	capture qui reg dv announce##c.time_a [pweight=weight]
	restore

	*if not missing observations
	if missing("`e(N)'")==0{
		qui lincom 1.announce 
		matrix results[`t',1]=r(estimate)
		scalar p = 2*ttail(r(df),abs(r(estimate)/r(se)))
		matrix results[`t',2]=p
		matrix results[`t',3]= r(estimate) - invt(r(df),.975)*r(se)
		matrix results[`t',4]= r(estimate) + invt(r(df),.975)*r(se)
		matrix results[`t',5]=`t'
		matrix results[`t',6] = bw
		matrix results[`t',7] = `max_poss'
	} 
	else if missing("`e(N)'")==1 {
		disp in red "No Observations"
	}
	
}


capture drop results*
svmat results	  

drop x

*add noise to the max possible sentence so that all the results are not stacked against each other
gen x = log(results7+rnormal(0,50))
twoway (rspike  results3 results4 x,lcolor(maroon) horizontal) ///
		(scatter x results1, mcolor(gs6) msymbol(circle_hollow) msize(0.5)), ///
       yscale() name(gr_pvalues_unadj, replace) ytitle("Charge Severity (log days)", size(medium)) ///
	   xline(0, lcolor(black) lpattern(dash)) ///	   
	   xtitle("Estimate and 95% CI", size(medium)) legend(off) ///
	   ylabel(6.8(.5)9, angle(0) nogrid) xlab(,nogrid) ///
	   title("Charge-FE RD Estimates and 95% Confidence Intervals, Unadjusted",size(medium))  graphregion(color(white)) 
graph export "charge_balance_count.pdf", replace

/********************************************************************************
*** Figure B.2. Tests for Bandwidth Artifacts - Petition Only 				  ***
********************************************************************************/
matrix bw_results_unadj = J(81,4,.)
matrix bw_results_adj = J(81,4,.)
forval bw = 10(1)90 {
	rdrobust sentence_normed_truncated time_a if abs(time_a)<180 , h(`bw' `bw')  b(`bw' `bw') ///
			  kernel(tri) vce(cluster cluster_county_statute) all
	matrix Vl = e(V_rb_l)
	matrix Vr = e(V_rb_r)
	local se = sqrt(Vl[1,1]+Vr[1,1])
	matrix bw_results_unadj[`bw'-9,1]=`bw'
	matrix bw_results_unadj[`bw'-9,2]=e(tau_bc) 
	matrix bw_results_unadj[`bw'-9,3]=e(tau_bc) - 1.959964*`se'
	matrix bw_results_unadj[`bw'-9,4]=e(tau_bc) + 1.959964*`se'

	rdrobust sentence_normed_truncated time_a if abs(time_a)<180 , h(`bw' `bw')  b(`bw' `bw') ///
			  covs(numberofcounts s_* j_*) kernel(tri) vce(cluster cluster_county_statute) all
	matrix Vl = e(V_rb_l)
	matrix Vr = e(V_rb_r)
	local se = sqrt(Vl[1,1]+Vr[1,1])
	matrix bw_results_adj[`bw'-9,1]=`bw'
	matrix bw_results_adj[`bw'-9,2]=e(tau_bc) 
	matrix bw_results_adj[`bw'-9,3]=e(tau_bc) - 1.959964*`se'
	matrix bw_results_adj[`bw'-9,4]=e(tau_bc) + 1.959964*`se'
}	
capture drop bw_results*
svmat bw_results_unadj
svmat bw_results_adj
local ar = 0.67
twoway (line bw_results_unadj3 bw_results_unadj1, lcolor(black) lpattern(dot)) ///
	   (line bw_results_unadj4 bw_results_unadj1, lcolor(black) lpattern(dot)) ///
	   (line bw_results_unadj2 bw_results_unadj1, lcolor(black) lpattern(solid)), ///
	   legend(off) name(gr_band_unadj, replace) xtitle("Bandwidth") ytitle("Estimate and 95% CI") /// 
	   xlab(10(10)90, nogrid) ylab(-0.3(0.1)0.5, nogrid format(%3.1f)) yline(0,lpattern(solid)) xline(45.1,lpattern(solid)) ///
	   aspectratio(`ar') title("Unadjusted Estimates, Varying Bandwidth") 

twoway (line bw_results_adj3 bw_results_adj1, lcolor(black) lpattern(dot)) ///
	   (line bw_results_adj4 bw_results_adj1, lcolor(black) lpattern(dot)) ///
	   (line bw_results_adj2 bw_results_adj1, lcolor(black) lpattern(solid)), ///
	   legend(off) name(gr_band_adj, replace) xtitle("Bandwidth") ytitle("Estimate and 95% CI") /// 
	   xlab(10(10)90, nogrid) ylab(-0.3(0.1)0.5, nogrid format(%3.1f)) yline(0,lpattern(solid)) xline(44.6,lpattern(solid)) ///
	   aspectratio(`ar') title("Adjusted Estimates, Varying Bandwidth") 

gr combine gr_band_unadj gr_band_adj, imargin(0 0 0 0) graphregion(margin(l=22 r=22)) xsize(7) ysize(3.5)
graph export "bandwidth_tests.pdf", replace		


  
/********************************************************************************
*** Figure B.3. Replication using Washington State Data 					  ***
********************************************************************************/
*The following section uses proprietary data initially collected from the Washington State Department of Corrections (DOC). See the "Source Data\Additional Data\WA Trial Courts" folder for instructions on generating these data from the DOC records
  
clear
use "recall_data_washington.dta", 

**generate calendar
gen time_a = edate -  mdy(6,6,2016)
gen announce = time_a>=0
egen cluster_judge_statute =group(judge Statute_RCW)
set matsize 10000
encode Statute_RCW, gen(charge_code)
encode judge, gen(judge_code)
recode TOTALCONFINEMENTDAYS .=0
local vrange = "0.2(0.1)0.6"
local fitmodel = "qfit"
local fitmodel2 = "lpoly"
local ar = 0.67

*** a. Unadjusted  

preserve
rdbwselect sentence_normed_truncated time_a if abs(time_a)<180, kernel(tri) vce(cluster cluster_judge_statute)	
scalar bw = e(h_mserd)
keep if abs(time_a)<=bw
fastxtile n_all = time_a, nq(50)
gen n = .
replace n = n_all
collapse (mean) sentence_normed_truncated (median) time_a, by(n)
gen weight = max(0,(bw-abs(time_a))/bw)
twoway 	   (`fitmodel2' sentence_normed_truncated time_a if time_a <= 0, lcolor(gs8) lpattern(solid)) ///
	       (`fitmodel2' sentence_normed_truncated time_a if time_a >=0, lcolor(gs8) lpattern(solid)) ///
		   (`fitmodel' sentence_normed_truncated time_a if time_a <= 0 [aweight=weight], lcolor(maroon) lpattern(dash)) ///
	       (`fitmodel' sentence_normed_truncated time_a if time_a >=0 [aweight=weight], lcolor(maroon) lpattern(dash)) ///
		   (scatter sentence_normed_truncated time_a, mcolor(black) msymbol(circle_hollow)), ///
			name(graph_announce_unadjusted_wa, replace) ///
			title("Petition Announcement, Unadjusted")  xtitle("Date") ///
			xlabel( 0 "Jun 6" -45 "Apr 22" 45 "Jul 21",angle(45) nogrid) ///
			ytitle("Normalized Sentence") ylabel(`vrange', format(%3.1f) nogrid) ///
			xline(0, lcolor(black) lpattern(dash)) legend(off) scheme(plotplain) aspectratio(`ar') ///
			graphregion(color(white)) bgcolor(white)
graph export "WA_announce_normal_no_control.pdf", replace
restore

*** b. Adjusted    
         										      
preserve	
rdbwselect sentence_normed_truncated time_a if abs(time_a)<180, kernel(tri) vce(cluster cluster_judge_statute)	
scalar bw = e(h_mserd)
qui reg sentence_normed_truncated i.judge_c i.charge_code  if abs(time_a)<=bw
qui predict ry if e(sample), resid
qui summ sentence_normed_truncated if e(sample), meanonly
qui replace ry = r(mean) + ry
qui reg time_a i.judge_c i.charge_code if abs(time_a)<=bw
qui predict rx if e(sample), resid
qui summ time_a if e(sample), meanonly
qui replace rx = r(mean)+rx 
fastxtile n_all = rx, nq(50)
gen n = .
replace n = n_all
collapse (mean) ry (median) rx time_a, by(n)
gen weight = max(0,(bw-abs(rx))/bw)
sum ry rx
twoway 	   (`fitmodel2' ry rx if rx <= 0, lcolor(gs8) lpattern(solid)) ///
	       (`fitmodel2' ry rx if rx  >=0 , lcolor(gs8) lpattern(solid)) ///
	       (`fitmodel' ry rx if rx <= 0 [aweight=weight], lcolor(maroon) lpattern(dash)) ///
	       (`fitmodel' ry rx if rx  >=0 [aweight=weight], lcolor(maroon) lpattern(dash)) ///
		   (scatter ry rx, mcolor(black) msymbol(circle_hollow)) ///
		    ,  name(graph_announce_adjusted_wa, replace) ///
			title("Petition Announcement, Adjusted")  xtitle("Date") ///
			xlabel( 0 "Jun 6" -45 "Apr 22" 45 "Jul 21",angle(45) nogrid) ///
			ytitle("Normalized Sentence") ylabel(`vrange', format(%3.1f) nogrid) ///
			xline(0, lcolor(black) lpattern(dash)) legend(off) scheme(plotplain) aspectratio(`ar') ///
			graphregion(color(white)) bgcolor(white)
			graph export "WA_announce_normal_control.pdf", replace
restore

gr combine graph_announce_unadjusted_wa graph_announce_adjusted_wa, ///
		   imargin(0 0 0 0) graphregion(margin(l=22 r=22)) name(all_rd_plots, replace)

graph export "combined_rd_plots_wa.pdf", replace
	  

/********************************************************************************
*** Table C.1. Hetereogenous Effects by Leniency - Petition Only 			  ***
********************************************************************************/
use "recall_data_main.dta", clear
*residualize normalized sentence before the petition announcement
reghdfe sentence_normed_truncated  numberofcounts if time_a<0, absorb(statute_code judge_code, savefe)

*generate indicator whether judge is above median severity
sum  __hdfe2__, d //the coefficients associated with the second fixed effect i.e. judges
gen x = __hdfe2__>`r(p50)'
replace x = . if mi(__hdfe2__)
bys judge: egen above_median = max(x) 

*** a. Above Median Leniency 

matrix tablec1 = J(6,6,.)
rdrobust sentence_normed_truncated time_a if abs(time_a)<180 & above_median==1, kernel(tri) vce(cluster cluster_judge_statute)  all
matrix b_temp = e(beta_p_l)
matrix v_temp = e(V_rb_l)
matrix tablec1[1,1]=  e(tau_bc)
matrix tablec1[2,1]= e(se_tau_rb)
matrix tablec1[3,1]=b_temp[1,1]
matrix tablec1[4,1]=sqrt(v_temp[1,1])		
matrix tablec1[5,1]=e(b_l)
matrix tablec1[6,1]=e(N_h_l)+e(N_h_r)
rdrobust sentence_normed_truncated time_a if abs(time_a)<180 & above_median==1,  kernel(tri) vce(cluster cluster_judge_statute) 	covs(numberofcounts s_* j_*) all
matrix b_temp = e(beta_p_l)
matrix v_temp = e(V_rb_l)
matrix tablec1[1,2]=  e(tau_bc)
matrix tablec1[2,2]= e(se_tau_rb)
matrix tablec1[3,2]=b_temp[1,1]
matrix tablec1[4,2]=sqrt(v_temp[1,1])		
matrix tablec1[5,2]=e(b_l)
matrix tablec1[6,2]=e(N_h_l)+e(N_h_r)

*** b. Below Median Leniency                                             

rdrobust sentence_normed_truncated time_a if abs(time_a)<180 & above_median==0, kernel(tri) vce(cluster cluster_judge_statute)  all
matrix b_temp = e(beta_p_l)
matrix v_temp = e(V_rb_l)
matrix tablec1[1,3]=  e(tau_bc)
matrix tablec1[2,3]= e(se_tau_rb)
matrix tablec1[3,3]=b_temp[1,1]
matrix tablec1[4,3]=sqrt(v_temp[1,1])		
matrix tablec1[5,3]=e(b_l)
matrix tablec1[6,3]=e(N_h_l)+e(N_h_r)
rdrobust sentence_normed_truncated time_a if abs(time_a)<180 & above_median==0,  kernel(tri) vce(cluster cluster_judge_statute) covs(numberofcounts s_* j_*) all
matrix b_temp = e(beta_p_l)
matrix v_temp = e(V_rb_l)
matrix tablec1[1,4]=  e(tau_bc)
matrix tablec1[2,4]= e(se_tau_rb)
matrix tablec1[3,4]=b_temp[1,1]
matrix tablec1[4,4]=sqrt(v_temp[1,1])		
matrix tablec1[5,4]=e(b_l)
matrix tablec1[6,4]=e(N_h_l)+e(N_h_r)

/*This is Table C1*/
matrix list tablec1

/********************************************************************************
*** Table C.2. Effect on Appellate Judges									  ***
********************************************************************************/

use "recall_data_appellate.dta", clear
destring _all, replace
egen panel_code = group(panel)
qui tab panel_code, gen(j_)

egen appeal_division_code = group(appeal_division)
qui tab appeal_division, gen(d_)

capture drop time_a 
capture drop time_r
gen time_a = date(edate,"YMD")-date("06/06/2016","MDY")
gen time_r = date(edate,"YMD")-date("06/05/2018","MDY")

capture drop table*
matrix tablec2 = J(6,4,.)

*** a. Petition Announcement

rdrobust prodefendant time_a if abs(time_a)<180, kernel(tri) vce(cluster panel_code) all
matrix b_temp = e(beta_p_l)
matrix v_temp = e(V_rb_l)
matrix tablec2[1,1]=  e(tau_bc)
matrix tablec2[2,1]= e(se_tau_rb)
matrix tablec2[3,1]=b_temp[1,1]
matrix tablec2[4,1]=sqrt(v_temp[1,1])		
matrix tablec2[5,1]=e(b_l)
matrix tablec2[6,1]=e(N_h_l)+e(N_h_r)
rdrobust prodefendant time_a if abs(time_a)<180, kernel(tri) vce(cluster panel_code) covs(j_*) all
matrix b_temp = e(beta_p_l)
matrix v_temp = e(V_rb_l)
matrix tablec2[1,2]=  e(tau_bc)
matrix tablec2[2,2]= e(se_tau_rb)
matrix tablec2[3,2]=b_temp[1,1]
matrix tablec2[4,2]=sqrt(v_temp[1,1])		
matrix tablec2[5,2]=e(b_l)
matrix tablec2[6,2]=e(N_h_l)+e(N_h_r)

*** b. Recall Election

rdrobust prodefendant time_r if abs(time_r)<180, kernel(tri) vce(cluster panel_code) all
matrix b_temp = e(beta_p_l)
matrix v_temp = e(V_rb_l)
matrix tablec2[1,3]=  e(tau_bc)
matrix tablec2[2,3]= e(se_tau_rb)
matrix tablec2[3,3]=b_temp[1,1]
matrix tablec2[4,3]=sqrt(v_temp[1,1])		
matrix tablec2[5,3]=e(b_l)
matrix tablec2[6,3]=e(N_h_l)+e(N_h_r)
rdrobust prodefendant time_r if abs(time_r)<180, kernel(tri) vce(cluster panel_code) covs(j_*) all
matrix b_temp = e(beta_p_l)
matrix v_temp = e(V_rb_l)
matrix tablec2[1,4]=  e(tau_bc)
matrix tablec2[2,4]= e(se_tau_rb)
matrix tablec2[3,4]=b_temp[1,1]
matrix tablec2[4,4]=sqrt(v_temp[1,1])		
matrix tablec2[5,4]=e(b_l)
matrix tablec2[6,4]=e(N_h_l)+e(N_h_r)

/*This is Table C2*/
matrix list tablec2

/********************************************************************************
*** Figure C.3.	Effect on Sentencing of Signature Certification	Date		  ***
********************************************************************************/
*see lines 26 to 195